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1 Introduction 

The complementarity problem has become one of the most well-established disci¬ 
plines within mathematical programming [10], in the last three decades. It is not 
snrprising that the complementarity problem has received much attention of re¬ 
searchers, due to its widespread applications in the fields of engineering, economics 
and sciences. In the literature, many theoretical results and efficient numerical 
methods were developed, we refer the reader to [11] for an exhaustive survey on 
complementarity problems. 

The eigenvalue complementarity problem (EiCP) not only is a special type of 
complementarity problems, but also extends the classical eigenvalue problem which 
can be traced back to more than 150 years (see [12,30]). EiCP first appeared in 
the study of static equilibrium states of mechanical systems with unilateral friction 
[8], and has been widely studied [1,9,14,15,16] in the last decade. Mathematically 
speaking, for two given square matrices A,B£ EiCP refers to the task of 

finding a scalar A G M and a vector x G IR"'\{0} such that 

0 < a; _L to := {XB — A)x > 0. 

EiCPs are closely related to a class of differential inclusions with nonconvex pro¬ 
cesses defied by linear complementarity conditions, which serve as models for many 
dynamical systems. Given a linear mapping A : R" —^ R", consider a dynamic sys¬ 
tem of the form: 

f u(t) > 0, 

< u{t) — Au(t) > 0, (1-1) 

[ {u{t),ii(t) — Au{t)) = 0. 

It is obvious that (1.1) is equal to u{t) G F{u(t)), where the process F : R" —>• R" 
is given by 

Gr{F) := { (i, 2/) G R” x R” j x > 0,y — Ax > 0, {x, y — Ax) =0} 

and is nonconvex. As noticed already by Rockafellar [26], the change of variable 
u{t) = e^^v{t) leads to an equivalent system 

Xv{t) + v{t) G F{v{t)). 

This transformation efficiently utilizes the positive homogeneity of F. Therefore, 
if the pair (A, a;) satisfies Xx G F{x), then the trajectory t i->- e^*x is a solution of 
dynamic system (1.1). Moreover, if such a trajectory is nonconstant, then x must 
be a nonzero vector, which further implies that (A, x) is a solution of EiCP with 
B := I (i.e., B is the identity matrix). The reader is referred to [8,27] for more 
details. 



On the cone eigenvalue complementarity problem for higher-order tensors 


3 


When B is symmetric positive definite and A is symmetric, EiCP is symmetric. 
In this case, it is well analyzed in [25] that EiCP is equivalent to finding a sta¬ 
tionary point of a generalized Rayleigh quotient on a simplex. Generally speaking, 
the resulting equivalent optimization formulation is NP-complement [6,25] and 
very difficult to be solved efficiently, and in particular when the dimension of the 
problem is large. 

In the current numerical analysis literature, considerable interest has arisen in 
extending concepts that are familiar from linear algebra to the setting of multilin¬ 
ear algebra. As a natural extension of the concept of matrices, a tensor, denoted 
by A, is a multidimensional array, and its order is the number of dimensions. 
Let m and n be positive integers. We call A = ( 0 ^^...^^), where € R for 

1 < ii,...,im < n, a real m-th order n-dimensional square tensor. The tensor A is 
further called symmetric if its entries are invariant under any permutation of their 
indices. The eigenvalues and eigenvectors of such square tensor were introduced 
by Qi [20], and were introduced independently by Lim [18]. 

For a vector x = (xi,... ,Xn)^ G C", Ax'^~^ is an n-vector with its i-th com¬ 
ponent defined by 

n 

{Ax""~^), = ^ au 2 --i^Xi^---Xi^, for i = 1,2,... ,n, 

^2v,*m, = l 

and Ax^ is a homogeneous polynomial defined by 

n 

Ax ^ ^ ^Zi22---2r7i ^^2 ' ' ' ^*T7i ■ 

For given tenors A and B with same structure, we say that {A^B) is an identical 
singular pair, if 

|x G C"\{0} : Ax^-^ = = o| 0. 

Definition 1.1 ([5]) Let A and B be two m-th order n-dimensional tensors on R. 
Assume that {A,B) is not an identical singular pair. We say (A, x) G C x (C"\{0}) 
is an eigenvalue-eigenvector of {A,B), if the n-system of equations: 

(A-AR)x"""^ = 0, (1.2) 

that is, 

n 

^ ^ ^ - 1 , 2 ,...,?!, 

possesses a nonzero solution. Here, A is called a ^-eigenvalue of A^ and x is called 
a ^-eigenvector of A. 
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With the above definition, the classical higher order tensor generalized eigenvalue 
problem (TGEiP) is to hnd a pair of (A,®) satisfying (1.2). It is obvious that if 
B = I, the unit tensor I = (5^^...^^), where is the Kronecker symbol 

jl, if ii = • • • = im^ 

I 0, otherwise, 

then the resulting ^-eigenvalues reduce to the typical eigenvalues, and the real B- 
eigenvalues with real eigenvectors are the iil-eigenvalues, in the terminology of [20, 
22]. In the literature, we have witnessed that tensors and eigenvalues/eigenvectors 
of tensors have fruitful applications in various fields such as magnetic resonance 
imaging [3,24], higher-order Markov chains [19] and best-rank one approximation 
in date analysis [23], whereby many nice properties such as the Perron-Frobenius 
theorem for eigenvalues/eigenvectors of nonnegative square tensor have been well 
established, see, e.g., [4,31]. 

In this paper, we consider the tensor generalized eigenvalue complementarity prob¬ 
lem (TGEiCP), which can be mathematically characterized as finding a nonzero 
vector X G R” and a scalar A £ R with property 

x£K, XBx^~^ - Ax'^~^ € K*, {x,XBx"'~^ - Ax'^~^) = 0, (1.3) 

where A and B are two given m-th order n-dimensional higher tensors, is a closed 
and convex cone in R", and K* is the positive dual cone of /f, i.e., K* := {w € 
R” : (w,fc) > 0,V fc € ff}. As EiGPs closely relate to differential inclusions with 
processes defined by linear complementarity conditions, TGEiCPs are also closely 
related to a class of differential inclusions with nonconvex processes H defined by 

Gt{H) := {{x,y) £ R"xR" j x € K,By™'~^ —Ax"^~^ £ K*, {x,By'^~^—Ax'^~^) = 0}. 

The scalar A and the nonzero vector x satisfying system (1.3) are respectively 
called a A-eigenvalue of (A, B) and an associated A-eigenvector. In this situation, 
(A, x) is also called a A-eigenpair of {A,B). The set of all eigenvalues is called the 
A-spectrum of {A,B), and it is dehned by 

aK{A,B) := {X£R : 3x £ R^ViO}, A 9 x _L AZ3x™"^ - £ A*}. 

Throughout this paper one assumes that A n (—A) = {0} and Sx™ 0 for any 
X £ A\{0}. If A = {x £ R" : X > 0}, then (1.3) reduces to 

X > 0, XBx™'~^ — Ax^~^ >Q, {x,XBx™'~^ — Ax^~^) = 0, (1-4) 

which is a specialization of TGEiP. The scalar A and the nonzero vector x satisfy¬ 
ing system (1.4) are called a Pareto-eigenvalue of {A,B) and an associated Pareto- 
eigenvector, respectively. The set of all Pareto-eigenvalues, defined by a{A,B), is 
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called the Pareto-spectmm of If in addition m = 2, the problem under 

consideration immediately reduces to the classical EiCP. If ® € int(-K') (respec¬ 
tively, ® € {a; € R" : x > 0}), then A is called a strict Jf-eigenvalue (respec¬ 
tively, Pareto-eigenvalue) of In particular, if S = I, then the K (Pareto)- 

eigenvalue/eigenvector of {A,B) is called the K (Pareto)-eigenvalue/eigenvector of 
A^ and the K (Pareto)-spectrum of [A, B) is called the K (Pareto)-spectrum of A. 

The main contributions of this paper are four folds. As we have mentioned in 
above, TGEiCP is an essential extension of EiCP. Accordingly, a natural question 
is that whether TGEiCP has solutions like EiCP. In this paper, we first give an 
affirmative answer to this question, thereby discussing the existence of the solution 
of TGEiCP (1.3) under some conditions. Note that TGEiCP is also a special case 
of complementarity problem, and it is well documented in [10] that one of the 
most popular avenues to solve complementarity problems is reformulating them as 
optimization problems. Hence, we here also introduce two equivalent optimization 
reformulations of TGEiCP, which further facilitates the analysis of upper bound 
of cone eigenvalues of tensor. With the existence of the solution of TGEiCP, ones 
may be further interested in a truth that how many eigenvalues exist. Therefore, 
the third objective of this paper is to establish theoretical results concerning the 
bounds of the number of eigenvalues of TGEiCP. Finally, we develop a projection 
algorithm to solve TGEiCP, which is an easily implementable algorithm as long 
as the convex cone K is simple enough in the sense that the projection onto K has 
explicit representation. As an illustration of our theoretical results, we implement 
our proposed projection algorithm to solve some synthetic examples and report 
the corresponding computational results. 

The structure of this paper is as follows. In Section 2, the existence of solution 
for TGEiCP is discussed under some reasonable assumptions. Two optimization 
reformulations of TGEiCP are presented in Section 3, and the relationship of 
TGEiCP with the optimization of the Rayleigh quotient associated to tensors is 
established. Moreover, based upon a reformulated optimization model, an upper 
bound of cone eigenvalues of tensor is also established. In Section 4, some the¬ 
oretical results concerning the bounds of number of eigenvalues of TGEiCP are 
presented. To solve TGEiCP, we develop a so-called scaling-and-projection algorithm 
(SPA) and conduct some numerical simulations to support our results of this pa¬ 
per. Finally, we complete this paper with drawing some concluding remarks in 
Section 6. 

Notation. Let M" denote the real Euclidean space of column vectors of length n. 
Denote R" = {x € R” : a; > 0} and R"_|_ = {a; € R" : a; > 0}. Let A be a tensor of 
order m and dimension n, and J be a subset of the index set N := {1, 2,..., n}. We 
denote the principal sub-tensor of A by Aj, which is obtained by homogeneous 
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polynomial Ax^ for all x = (xi, X2, ■ ■ ■, Xn)^ with Xi = 0 for N\J. So, Aj is a 
tensor of order m and dimension | J|, where the symbol \J\ denotes the cardinality 
of J. For a vector x € K" and an integer r > 0, denote = (xi, X 2 ,. •., x^)^. 


2 Existence of the solution for TGEiCP 


This section deals with the existence of the solution for TGEiCP. Let if be a closed 
and convex pointed cone in M". Recall that a nonempty set S C M" generates a 
cone K and write K := cone{S') if K := {ts : s G S, t € R+}. If in addition S does 
not contain zero and for each k G if\{0}, there exists unique s € S and t G 1R+ 
such that k = ts, then we say that S' is a basis of K. Whenever S is a hnite set, 
cone(conv(S)) is called a polyhedral cone, where conv(S) stands for the convex hull 
of S. Let if be a closed convex cone equipped with a compact basis S. To study 
the existence of solution for TGEiCP, we first make the following assumption. 

Assumption 2.1 It holds that Sx™ ^ 0 for every vector x G S. 


Remark 2.1 It is easy to see that Assumption 2.1 holds if and only if one of the 
tensors B (or —B) is strictly if-positive, i.e., Bx^ > 0 (or —Rx™ > 0) for any 
X G if\{0}. In particular, when if = R", if R is a strictly copositive tensor (see 
[21,29]), then B satishes Assumption 2.1. It is easy to see that if B is nonnegative, 
i.e., S > 0, and there are no index subset J of A such that Bj is a zero tensor, 
then Sx'" > 0 for any x G R!jl\{0}, and hence, in this case. Assumption 2.1 holds. 


From (1.3), one knows that if (A,x) G R x (R"\{0}) is a if-eigenpair of {A,B), 
then necessarily 

Sx™ ’ 

provided Rx"* ^ 0. Consequently, by the second expression of (1.3), it holds that 


Ax'" 

Rx™ 


Bx^~ 


1 


-Ax" 


G if*. 


We now present the existence theorem of TGEiCP, which is a particular instance 
of Theorem 3.3 in [17]. However, for the sake of completeness, here we still present 
its proof. 


Theorem 2.1 Let K be a cone equipped with convex compact basis S. If Assumption 
2.1 holds, then TGEiCP (1.3) has at least one solution. 


Proof Define F : S x S' —>■ R by 


Ax'" 


{Bx”^-\y)- 


F{x,y) = {Ax'^ \t/) 


Rx'" 


( 2 . 1 ) 
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Since Bx^ 0 for any a; £ S', it is obvious that F{-,y) is lower-semicontinuous 
on S for any fixed y € S, and F{x, •) is concave on S for any fixed x € S. By the 
well-known Ky Fan inequality [2], there exists a vector x € S C -K'\{0} such that 

sup F{x, y) < sup F{y, y). (2.2) 

yes yes 

Consequently, since F{y, y) = 0 for any y £ S, by (2.2) it holds that F{x, y) < 0 for 
any y £ S. Let A = . Then, by (2.1), one knows that (\Bx"^~^ —Ax"^~^, y) > 0 

for any y £ S, which implies 


XBx^-^ - Ax"^-^ £ K*, 


(2.3) 


since for any y £ if it holds that y = ts for some t £ R-i- and s £ S. Moreover, it is 
easy to know that 

{x,XBx"^~^ - Ax^~^} =0, 

which means, together with (2.3) and the fact that x £ lsr\{0}, that (X,x) is a 
solution of (1.3). We obtain the desired result and complete the proof. □ 

From Theorem 2.1, we obtain the following corollary. 

Corollary 2.1 /f B is strictly copositive, then (1.4) has at least one solution. 

Proof Take S := {x £ R” | X]r=i = !}• i® clear that S' is a convex compact 
base of M" . By Theorem 2.1, it follows that the conclusion holds. The proof is 
completed. □ 


The following example shows that Assumption 2.1 is necessary to ensure the 
existence of the solution of TGEiCP. 


Example 2.1 Let m = 2. Consider the case where 


A = 



and 



It is easy to see that Assumption 2.1 does not hold for the above two matrices. 
Since det(AB — A) = —X^ — 11^0 for any A £ M, we claim that the system of 
linear equations [XB — A)x = 0 has only one unique solution 0 for any A £ R, 
which means that (A,r) £ R x R++ satisfying (1.4) does not exist. Moreover, we 
may check that {XB — A)x > 0 does not hold for any (A,a:) £ R x (R^\{0}) with 
X = (*1,0)^ or a; = (0,X2)^- Therefore, problem (1.4) has no solution. 
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3 Optimization reformulations of TGEiCP 

In this section, we study two optimization reformulations of (1.4). We begin with 
introducing a so-called generalized Rayleigh quotient related to tensors. For two 
given m-th order n dimensional tensors A and B, the related Rayleigh quotient is 
defined by 

Ax^ 

A(*) = (3.1) 

where Bx^ ^ 0. li m = 2, then X{x) defined by (3.1) reduces to one introduced in 
[25]. When A is symmetric and B is symmetric and strictly copositive, it is easy 
to see that the gradient of X{x) is 

(3-2) 

Notice that the expression (3.2) of the gradient of the Rayleigh quotient is only 
valid when A and B are both symmetric. Moreover, in this case, the stationary 
points of A(a:) correspond to solutions of (1.4). If either M or R is not symmetric, 
the above expression of VA(a:) is incorrect, and the relationship between stationary 
points and solutions of the TGEiCP with K = R" ceases to hold. 

The following lemma presents two fundamental properties of the generalized 
Rayleigh quotient A in (3.1), whose matrix version was proposed in [25]. Its proof 
is straightforward and skipped here. 

Lemma 3.1 For all x £ ]R"'\{0}, the following statements hold: 

(1) . X{tx) = X{x), Vr > 0; 

(2) . x^\/X{x) = 0. 

We first consider the following optimization problem 

p{A, B) := max { A(a;) j a: £ 5 } , (3.3) 

X 

where A(a;) is defined in (3.1), and the constraint set S is determined by 

5':=|a;£R+ : = l|, (3.4) 

which is called the standard simplex in M". 

We generalize the result of symmetric EiCP studied in [25] to TGEiCP as the 
following proposition. 

Proposition 3.1 Assume that the tensors A and B are symmetric and B is strictly 
copositive. Let x be a stationary point of (3.3). Then (A(®), x) is a solution of TGEiCP 
with if = R![l. 
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Proof Since a; is a stationary solution of (3.3), from the structure of S, there exist 
a G R" and d G R, such that 


—VA(2;) = a-\- ^e, 
a > 0,x > 0, 

X = 0 , 

e^a: = 1, 


(3.5) 


where e G R" is a vector of ones. By (3.5), we know —a;^VA(ai) = d, which implies, 
together with Lemma 3.1 (2), that j3 = 0. Consequently, from (3.2), the first two 
expressions of (3.5) and the fact that Bx"^ > 0, it holds that \{x)Bx'^~^ — Ax'^~^ > 
0. This means, together with the fact that s > 0 and x^{X{x)Bx^~^ —Ax'^~^) = 0, 
that (A(a:),x) is a solution of TGEiCP with K = R” . We complete the proof. □ 

In what follows, we denote 

A( 4 i^g = max {a : 3 a; G R" \{0} suet that (A, x) is a solution of (1.4)} 

for notational simplicity. Then, the following theorem characterizes the relation¬ 
ship between problem (3.3) and TGEiCP with K := R". 

Theorem 3.1 Let A and B be two m-th order n dimensional symmetric tensors. If B 
is strictly copositive, then = p{A,B). 

Proof It is obvious that the constrained set 17 of (3.3) is compact, and hence there 
exists a vector ® G 17 such that p{A, B) = \{x). It is clear that {e} U {e^ : i G I{x)} 
is linearly independent since x ^ 0, where I(x) = {i ^ N : x^ = 0}. Consequently, 
the hrst order optimality condition of (3.3) holds, which means that x is stationary 
point of (3.3). By Proposition 3.1, we know that (A(a;),a:) is a solution of TGEiCP 
with K = R". Hence, it holds that p{A,B) < A^^^g. 

Let (A, x) be a solution of TGEiCP with K := R", then A = Ax'^ jBx™'. Taking 
y = xjie^ x') implies that y G O. By Lemma 3.1 (1), we know A = Ay^ jBy^, which 
implies that A < p{A, B) from the definition of p{A, B). So, we have A^^g <p{A,B). 
Therefore, we obtain the desired result and complete the proof. □ 

We now study another optimization reformulation of TGEiCP with K := R" . 
We consider the following optimization problem 


7(AI, B) = max \Ax^ I x € U] , 

X ^ ^ 


(3.6) 


where E := {x G R" : Bx"^ = 1} is assumed to be compact. 
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Remark 3.1 li B is strictly copositive, then we claim that S is compact. Indeed, if 
S is not compact, then there exists a sequence C S such that ^ cio 

as fc —> oo. Taking := II clearly shows G R" and ||y^^^|| = 1. 

Without loss of generality, we may assume that there exists a vector y € R" 
satisfying |ii/|j = 1, such that —>■ j/ as fc —>■ oo. On the other hand, we have 

_ i/||a;(^) 11™^ which implies By™' = 0. It contradicts to the fact that 
By™ > 0, since y G R" \{0}. 

For TGEiCP with K := R" and (3.6), we have the following theorem which 
can be proved by a similar way to that used in [28] . 

Theorem 3.2 Let A and B be two m-th order n dimensional symmetric tensors. If 
Bx™ > 0 for any x G R" \{0}, then = 7{.4, B). 

By Theorems 3.1 and 3.2, it follows that solving the largest Pareto eigenvalue 
of TGEiCP is an NP-hard problem in general, i.e., there are no polynomial-time 
algorithm for solving the largest Pareto eigenvalue of TGEiCP. In the rest of this 
section, based upon Theorem 3.2, we further study the bound of Pareto eigenvalue 
of TGEiCP with B:=l and A" := R!]:. 

We denote by 12* the solution set of (1.4) with B := I and let 

= max •[ |A| : 3 a; G R+\{0} suet that (A, a:) G 12* } . 

Theorem 3.3 Suppose B : = T. It holds that 

|A|]4^^ < min|n^^|!7l|lF, a - , 

tu/iere a := max{ ^ \ < iiAi, ■ ■ ■ Am < n]. 


Proof Let (A, a;) be an arbitrary solution of (1.4) with B := X. Then it holds that 

^ ~ sr-n ’ 

Z^i=l ■''i 

which implies 

|7la:"^| ^ MIIfIIx^IIf 

l^i = \ '^i Z^2 = l *^2 

where := {xi^Xi^ • • • ^ which is an m-th order n-dimensional 

tensor. Since 

n / n 

II m\\2 \ ^ / \2 / \ ^ 2 \ ^ m—2 / \ m \ 

Ik ||f= 2^ ^ ) ’ 

2l,Z2,...,2m = l \2 = 1 / \z=l / 


|A| < n 2 MIIf- 


we obtain 
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On the other hand, we have 



m 


< a ■ n 


m— 1 


Hence we know 


|A| < min | 



Mil 


p,a ■ n 


m — 1 


} 


By the arbitrariness of A, we obtain the desired result and complete the proof. □ 

For the case where B is strict copositive bnt B hy & similar way, we may 
obtain 



where N^i^{B) = min{Bx”^ : x G R", = 1} > 0 by Theorem 5 in [21]. 


Notice that the computation of Nniin(H) is also NP-hard itself. 

4 Bounds for the number of Pareto eigenvalues 

In this section, we study the estimation of the nnmbers of Pareto-eigenvalue of 
{A,B)^ where A and B are two given m-th order n-dimensional tensors. We begin 
this section with some basic concepts and properties of eigenvalue/eigenvector of 
tensors. 

It is well known that, on the left-hand side of (1.2), {A—\B)x'^~^ is indeed a set 
of n homogeneous polynomials with n variables, denoted by {Pi{x) | 1 < i < n}, 
of degree (m — 1). In the complex field, to stndy the solution set of a system of 
n homogeneous polynomials (Pi,..., Pn), in n variables, the idea of the resultant 
Res (Pi,..., Pn) is well defined and introduced in algebraic geometry literatnre, 
we refer to the recent monograph [7] for more details. Applying it to our current 
problem, Res(Pi,..., Pn) has the following properties. 

Proposition 4.1 We have the following results: 

(1) . Res(Pi,... ,Pn) = 0, if and only if there exists {X,x) G C x (C"\{0}) such that 
satisfies (1.2). 

(2) . The degree of A in Res(Pi ,... ,Pn) is at most n(m — 1)"“^. 

For the considered TGEiCP with K = R" , we present the following proposition 
which fully characterizes the Pareto-spectrum of TGEiCP. 

Proposition 4.2 Let A and B be two m-th order n-dimensional tensors. A real num¬ 
ber \ is Pareto-eigenvalue of{A,B), if and only if there exists a nonempty subset J G N 
and a vector w G such that 


A m— 1 \ r? Tn— 1 

Ajw =XBjw 


(4.1) 
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and 


- au^...i^)wi^ ■ ■ ■ > 0, for every i £ N\J. (4.2) 


In such a case, the vector x € K” defined by 


Wi, i e J, 

0, i e N\J 


is a Pareto-eigenvector of (.4, B), associated to the real number \. 

Proof It can be proved by a similar way to that used in [28] and we skip it here. □ 

Remark f.l It is obvious that, in the case where B := T, (4.1) and (4.2) turn out 
to be 


\ m—1 \ fm—ll 

AjW = \w^ ^ 


(4.3) 


and 



(4.4) 


■ ■ ■ Wi,^ < 0, for every i £ N\J, 


respectively. The corresponding conclusions of Pareto-eigenvalues of A were stud¬ 
ied in [28]. 

By Proposition 4.2, if A is Pareto-eigenvalue of {A,B), then there exists a 
nonempty subset J C N such that A is a strict Pareto-eigenvalue of {Aj,Bj). 
Motivated by the works on estimating the cardinality of the Pareto-spectrum of 
matrices [27], we now state and prove the main results in this section. 

Theorem 4.1 Let A and B be two given m-th order n-dimensional tensors. Assume 
that {A,B) is not an identical singular pair. Then there are at most Sm,n ’.= 
Pareto-eigenvalues of{A,B). 

Proof It is obvious that, for every fc = 0, l,...,n — 1, there are („";!.) corresponding 
principal sub-tensors pair of order m dimension n — k. Moreover, by Proposition 
4.1, we know that every principal sub-tensors pair of order m dimension n — k can 
have at most (n — fc)(m — i)"--*:-i strict Pareto-eigenvalues. By Proposition 4.2, 
in this way one obtains the upper bound 



Hence proved. 


□ 
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Now we extend the above result to a more general case where A" is a polyhedral 
convex cone. A closed convex cone K in R" is said to be finitely generated if there 
is a linear independent collection {ci, C 2 , • • •, Cp} of vectors in R" such that 

AT = cone{ci, C 2 ,..., Cp} = QjCj : q = (ai, 02 ,..., G R^ | . (4.5) 

It is clear that K = {C^a \ a € R(J.}, where C = [ci,C 2 ,... ,Cp]^. Moreover, it is 
easy to see that the dual cone of K, denoted by A'*, K* = {w G R" | Cw > 0}. 

Theorem 4.2 Let A and B be two given m-th order n dimensional tensors. If the 
closed convex cone K admits representation (4.5), then {A,B) has at most 5m,p '■ = 
pm^~^ K-eigenvalues. 


Proof We hrst prove that problem (1.3) with K dehned by (4.5) is equivalent to 
finding a vector a G Rf'\{0} and A G R with property 

d>0, ADa"*"^> 0, (d,ADQ™"^= 0, (4.6) 


where P and Q are two m-th order p-dimensional tensors, whose elements are 
denoted by 

71 


and 

71 

jl,j2,...,jm = l 

for 1 < ii, 12 ,...,im < p, respectively. 

Let (x,A) G (R"\{0}) X R be a solution of (1.3) with K dehned by (4.5). Since 
a: G AT, by the dehnition of K, there exists a nonzero vector d G R(|. such that 
X = C^d. Consequently, from ABx”^~^ — Ax'"~^ G K* and the expression of K*, 
it holds that C{XBx^~^ — Ax"^~^) > 0, which implies 


C(AZ3(C'^d)'""^ - M(C'^d)”""^) > 0. (4.7) 


By the dehnitions of V and Q, we know that (4.7) can be equivalently written as 

AUd^-^-gd™-! >0. 


Moreover, it is easy to verify that (d, ADd"*“^ — C/d"*“^) = 0. Conversely, if (d, A) G 
(RP\{0}) X R satishes (4.6), then we can prove that (x. A) with x = a satishes 
(1.3) by a similar way. 

Consequently, by applying Theorem 4.1 to the problem (4.6), we know that 
{A,B) has at most 5m,p = pm‘‘~^ A'-eigenvalues. The proof is completed. □ 
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The above theorem shows that ax{A,B) has finitely many elements in case 
where 7^ is a polyhedral convex cone. However, in the nonpolyhedral case the 
situation can be even worse. For instance, lusem and Seeger [13] successfully con¬ 
structed a symmetric matrix A (i.e., 2-th order n dimensional tensor) and a non¬ 
polyhedral convex cone K such that ax {A, In) behaves like the Cantor ternary 
set, i.e., it is uncountable and totally disconnected. 

In the rest of this section, we discuss the case where B := I. We first present 
the following lemmas. 

Lemma 4.1 Let A be anm-th order n-dimensional nonnegative tensor, i.e., > 

0 for 1 < ii,... ,im ^ n. If A has two eigenvectors in then, the corresponding 

eigenvalues are equal. 

Proof Let Ai and A 2 be two Pareto-eigenvalues of A, and x,y £ K" + the corre¬ 
sponding associated Pareto-eigenvectors, which means 

= Aix[“-^1 and Ay^-^ = 

Since A is nonnegative tensor, we know that Ai, A 2 are nonnegative. Without loss 
of generality, assume Ai > A 2 . If Ai = 0, then A 2 = 0. Now we assume Ai > 0. 
Denote 

to = min{t >0 : ty — x £ R+j, (4-8) 

which must exist since y £ R" + . It is obvious that toy — x £ R", which implies 
that toyi > Xi for all i. Consequently, since > 0 for 1 < i\,... ,im < n, by 

the definitions of Ax"^~^ and A(toy)"‘~^, one knows that 

- Aia;["^-il = A(toy)”^-^ - Ax^-^ £ R!^, 

which implies 

to(A2/Ai)^^i/ - X £ R”. 

By (4.8), we know that to < to(A 2 /Ai) , which implies Ai < A 2 . Therefore, we 
obtain Ai = A 2 and complete the proof. □ 

Let A be an m-th order n-dimensional tensor, we say that .4 is a .Z-tensor, if 
all off-diagonal entries of A are nonpositive. 

Lemma 4.2 Let A be an m-th order n-dimensional tensor satisfying any of the fol¬ 
lowing conditions: (i) —A is a Z-tensor; (ii) A is a Z-tensor. Then, A admits at most 
one strict eigenvalue. 
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Proof We first consider case (i). Let Ai, A 2 £ R be two strict eigenvalues of A, i.e., 
there are vectors x,y ^ 1 ^++ such that and 


Hence, 


{A + y.l)x^ ^ = {\i+y)x^'^ and {A + yLT)y"^ ^ = (A 2 -|- 


where y is any real number. Since —A is a ^-tensor, A + yPL is nonnegative for y 
sufficiently large. By Lemma 4.1, we obtain the equality Xi + y = \2 + n, which 
implies the desired conclusion. 


In case (ii), the conclusion can be proved in a similar way. 


□ 


Proposition 4.3 Let A be an m-th order n-dimensional tensor satisfying any of the 
following conditions: (i) —A is a Z-tensor; (ii) A is a Z-tensor. Then, A can have at 
most pn := 2” — 1 Pareto eigenvalues. 

Proof We only consider case (i). The conclusion for case (ii) can be proved in a 
similar way. For every fc = 0,1,..., n — 1, there are {^-k) principal sub-tensors of 
order m dimension n — k. Since —A is a Z-tensor, it is clear that any principal 
sub-tensors of —A are also Z-tensor. Consequently, by Lemma 4.2, we know that, 
every principal sub-tensors can have at most one strict eigenvalues. Therefore, by 
Proposition 4.2, one gets the upper bound 



We obtain the desired result and complete the proof. 


□ 


It is easy to see that, if Al is a nonnegative tensor, then —A is a Z-tensor. Hence, 
by Proposition 4.3, we know that any m-th order n dimensional nonnegative tensor 
can have at most (2" — 1) Pareto eigenvalues. The following example shows that 
the bound pn is sharp within the second class mentioned in Proposition 4.3. This 
is what we call the exponential growth phenomenon. 

Example f.l Consider a 3-th order n-dimensional tensor A = (aiii 2 f 3 )i<ii,i 2 .i 3 <n 
with 0^^2213 = —and a > \fA. Given an arbitrary index set J = {Zi, I 2 , • • •, Ir} 
with 1 < l\ < I 2 < ■ ■ ■ < Ir < n, the principal sub-tensor Aj = {<^jih 33 )^< 3 i,j 2 ,j 3 <r 

1 I 1 I j Z Z 2 I -p ~v~ 

has jg = —a ^2 +‘ 23 . Take vector ( = (a^,a^,...,a^) . It is obvious that 




r 



r 



32,33 = ^ 


I2 J3 = l 




2 
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where Xj = — . This means that (4.3) holds. Since < 0 and 

^ > 0, one does not have to worry about the condition (4.4). By Remark 4.1, 
we know that \j is a Pareto-eigenvalue of A. Now we need to check that Aj^ ^ 
\j^ whenever Ji ^ J 2 . Take Ji,J 2 C {1,2, ...,n} with Ji ^ J 2 . Since J 1 AJ 2 = 
(di\J 2 ) U (J 2 \di) ^ 0, one can define k = max{fc € (1, 2,..., n}. A: e J 1 AJ 2 }. 
Without loss of generality, we assume that A: G J 2 , which implies A ^ Ji. In this 
case, we have 


leA jeJi j&J2,j<k jeJi,j<k-i 

3 

where b = . This implies that 


j&J2,j<k jeJi,j<k-l j = l 


fofe+i -2b’^ +b ^ b 
6-1 - 6-1 


> 0 , 


where the last inequality comes the fact 6 > 2 from the given condition that a > 
Therefore, we know that Aj^ ^ Aj^. Since there are 2" — 1 ways of choosing the 
index set J, there are as many elements in the Pareto spectrum of this special 
tensor A. 


Proposition 4.4 Suppose that there exists an index subset Jo ^ N with |Jo| = I 
such that > 0 for any i £ Jq and i 2 ,---,im € N\{i}. Then A has at most 

'1m,n '■= [n(m—1)+Z](m —l)*“^m”“*“^ Pareto-eigenvalues. In particular, if Jq =N, 
then A has at most pLm,n '■= n{m — 1)"“^ Pareto-eigenvalues. 

Proof Under the given condition, we only need to consider the principal sub-tensor 
Aj with Jo C J, which is due to the condition (4.2). Among the principal sub¬ 
tensors of order m dimension fc, there are (^Zj) of them with that property. This 
leads to the upper bound 

jL.n = E 

= {m- 1)' E s (s + 0("i - 1)*”^ 

= [n(m — 1) -|- Z](m — l)^“^m"“*“^. 

In particular, if Jq = N, we obtain immediately the desired result. The proof is 
completed. □ 

A similar type of argument leads to the following result: 

Proposition 4.5 Suppose that there exists an index set Jo Z A with |Jo| = I such 
that > 0 for any i £ Jq and i 2 , ■ ■ ■ ,im € JV\{i}. Moreover, suppose that —A 

is a Z-tensor. Then, A has at most ah := 2”“* Pareto-eigenvalues. 
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Proof This time one has to compute 


a 


l 

n 



n—l 

■i = E 


s=0 



• 1 = 2 ""'. 


We obtain the desired result and complete the proof. 


□ 


Theorems 4.1-4.2 and Propositions 4.3-4.5 extend the corresponding resnlts 
for bounds of Pareto eigenvalue of square matrix, which were studied in [27], to 
the case higher order tensors. In the square matrix case, i.e., m = 2, it is clear that 

^ Pn ^ 72,ra ^ ^2,71, 


which was presented in [27]. In the tensor case, i.e., m > 3, it is obvious that 
< pn and 7m,n < Sm,n for any 1 < Z < n. Moreover, it is not difficnlt to verify 
that, if Z = n then 7m,n = n{m — 1)"“^ > n2^~^ > pn', if 1 < Z < n — 1, then 
7m,n > [n{m — 1) + l](m — 1)”“^ > (2n + 1)2"“^ > pn- Therefore, it always holds 
that 

^ pn ^ 'ym,n ^ Sm,n 


for any 1 < Z < n. 


5 Numerical algorithm and simulations 

In this section, we first introduce an implementable algorithm for solving the 
TGEiCP. Then, we conduct some numerical results to verify the existence of the 
solntion of TGEiCP and the reliability of our proposed algorithm. 


5.1 Numerical algorithm 

It well known that the general nonlinear complementarity problem can also be 
transformed into a system of equations. Therefore, it is of course possible to ap¬ 
ply the semismooth and smoothing Newton methods to solve the problem under 
consideration in this paper. However, TGEiCP is more complicated than the clas¬ 
sical EiCP due to the high-dimensional structnre of tensor, thereby making snch 
second-order algorithms difficult to be implemented. Motivated by the recent work 
in [9] for solving matrix cone constrained eigenvalue problem, in this section, we 
extend the so-called scaling-and-projection algorithm (SPA), developed in [9], to 
solve (1.3) and follow the same name for TGEiCP. The corresponding algorithm 
can be described in Algorithm 1. Thronghout this section, we assume that B is 
strictly A-positive, i.e., Bx^ > 0 for any x € A\{0}. 
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Algorithm 1 A Scaling-and-Projection Algorithm (SPA). 

1: Take any starting point G /C\{0}, and define . 

2: for/c = 0,1, 2, ■ ■ ■ do 

3: One has a current point G iC\{0}. Compute 

= -d (5.1) 

4: If = 0, then stop. Otherwise, let := ||y^^^||, and compute 

(fc) 

= Uk and ^ (5.2) 

L J -y^(^{fc))m 

5: end for 


It is easy to verify that iterative scheme (5.1) always ensures = 0. 

As a consequence, € K* clearly means that {x^^\y^^'^) is a solution of problem 
(1.3). However, for the sake of convenience, we often use |j = 0 as the stopping 
condition in algorithmic framework instead of y^^^ G K*. 

As we have mentioned, our proposed algorithm is a straightforward extension of 
[9], we can easily get the following convergence theorem. For the sake of simplicity, 
we skip the corresponding proof of Algorithm 1, those who are interested in are 
referred to [9] for a similar proof. 

Theorem 5.1 Let the sequence he generated by Algorithm 1 and further satisfy 

g( 2 .(fe))"* = I Assume convergence of {x^^^} toward some limit that one denotes by 
X. Then, 

_ A _ 

lim Afe = A := ——, lim y^^^ = y := Ax"^~^ — \Bx^~^, (5.3) 

and (A, a;) is a solution of (1.3). 

Remark 5.1 As mentioned in [9], if K has a complicated structure, then computing 
in Algorithm 1 is not an easy task. However, there are many interesting cones 
for which the projection map admits an explicit and easily computable formula. 
This is true, for instance, for the Pareto cone, for the Loewner cone of positive 
semidefinite symmetric matrices, for the Lorentz cone and, more generally, for 
any revolution cone. Therefore Algorithm 1 is easily implemented as long as the 
projection onto K is easy enough to computed explicitly. 

Remark 5.2 The tensors A and B considered above are not necessarily symmetric. 
If A = R" and the tensors A and B are both symmetric, then the symmetric 
TGEiCP can be solved by computing a stationary point of the nonlinear program 
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(3.3). The constraint set of this program is the simplex S defined by (3.4). The 
special structure of this set S makes the computation of projections of vectors over 
S very easy. On the other hand, the objective fnnction of the reqnired nonlinear 
program has Hessian whose computation is quite involved. These features lead to 
our decision of investigating first order algorithms that are based on gradients and 
projections. This will be our investigation task in future. 


5.2 Numerical simulations 


We have theoretically discussed the existence of the solution of TGEiCP in Section 
2 and introduced an implementable projection method to solve the problem under 
consideration in Section 5.1. Thus, in this section, we aim at verifying that our 
theoretical results are true, in addition to demonstrating the reliability of the pro¬ 
posed algorithm. We implement Algorithm 1 by Matlab R2012b and conduct the 
numerical simulations on a Lenovo notebook with Inter(R) Core(TM) i5-2410M 
CPU 2.30GHz and 4GB RAM running on Windows 7 Home Premium operating 
system. 

In our experiments, we concentrate on three concrete TGEiCPs with symmetric 
structure and only list the details of tensors A and B in the ensuing examples. 


Example 5.1 We consider two 4-th order 2-dimensional symmetric tensors A and 
B, where the tensor A is specihed as 


A(:,:,l,l) = 

A(:,:,2,l) = 


0.8147 0.5164 
0.5164 0.9134 

0.4218 0.8540 


0.8540 0.9595 
and the tensor B is specihed as 

1.6324 1.1880 


H(:,:,l,l) = 

H(:,:,2,l) = 


1.1880 1.5469 

1.6557 1.4424 \ 
1.4424 1.9340 j 


A(:,:,l,2) = 

A(:,:,l,2) = 

H(:,:,l,2) = 

H(:,:,l,2) = 


0.4218 0.8540 \ 
0.8540 0.9595 


0.6787 0.7504 
0.7504 0.3922 


1.6557 1.4424 
1.4424 1.9340 

1.6555 1.4386 
1.4386 1.0318 


Example 5.2 This example considers two 4-th order 3-dimensional symmetric ten¬ 
sors A and B, where A and B are specihed as follows: 


A(:,:,l,l) 


/ 0.6229 0.2644 0.3567 \ 

0.2644 0.0475 0.7367 , A(:,:,l,2) 

V 0.3567 0.7367 0.1259/ 


/ 0.7563 0.5878 0.5406 \ 
0.5878 0.1379 0.0715 , 

V 0.5406 0.0715 0.3725/ 
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^(:,:,3,3) 

B{-.,-.,2,3) 

B{-.,:,3,2) 


/ 0.0657 
0.4918 

V 0.9312 

/ 0.7689 
0.3941 

V 0.6034 

/ 0.0657 
0.4918 

V 0.9312 

/ 0.7581 
0.7205 

V 0.9044 

/ 0.6730 
0.5351 

V 0.4473 

/ 0.6730 
0.5351 

V 0.4473 

/ 0.4632 
0.2043 

V 0.2823 

/ 0.4632 
0.2043 

V 0.2823 


0.4918 

0.7788 

0.9045 

0.3941 

0.3577 

0.3465 

0.4918 

0.7788 

0.9045 

0.7205 

0.0782 

0.7240 

0.5351 

0.2853 

0.3071 

0.5351 

0.2853 

0.3071 

0.2043 

0.7282 

0.7400 

0.2043 

0.7282 

0.7400 


0.9312 \ 
0.9045 
0.8711 / 

0.6034 \ 
0.3465 
0.4516 / 

0.9312 \ 
0.9045 
0.8711 / 

0.9044 \ 
0.7240 
0.3492 ) 

0.4473 \ 
0.3071 
0.9665 / 

0.4473 \ 
0.3071 
0.9665 / 

0.2823 \ 
0.7400 
0.9369 / 

0.2823 \ 
0.7400 
0.9369 / 


, A{:,-.,2,1) 

, A{:,-.,2,3) 

, A{:,-.,3,2) 

, B{-.,:AA) 

, 1,3) 

, B{-.,:,2,2) 

, 3,1) 

, B{:,:,3,3) 


0.7563 

0.5878 

0.5406 \ 

0.5878 

0.1379 

0.0715 

0.5406 

0.0715 

0.3725 / 

0.8077 

0.4910 

0.2953 \ 

0.4910 

0.5054 

0.5556 

0.2953 

0.5556 

0.9608 / 


/ 0.8077 0.4910 0.2953 
0.4910 0.5054 0.5556 

V 0.2953 0.5556 0.9608 

/ 0.6954 0.4018 0.1406 \ 
0.4018 0.9957 0.0483 , 

V 0.1406 0.0483 0.0988/ 

/ 0.7585 0.6433 0.2306 \ 
0.6433 0.8986 0.3427 , 

V 0.2306 0.3427 0.5390/ 

/ 0.3608 0.3914 0.5230 \ 
0.3914 0.6822 0.5516 , 

V 0.5230 0.5516 0.7091/ 

/ 0.7585 0.6433 0.2306 \ 
0.6433 0.8986 0.3427 , 

V 0.2306 0.3427 0.5390/ 

/ 0.8200 0.5914 0.4983 \ 
0.5914 0.0762 0.2854 . 

V 0.4983 0.2854 0.1266/ 



Example 5.3 This example also considers two 4-th order 3-dimensional symmetric 
tensors A and B, where A and B take their components as follows; 


- 4 (:,;, 1 , 1 ) 


-4(:,;,1,3) 


A{-.,:,2,2) 


0.4468 

0.4086 

0.5764 

0.4086 

0.8176 

0.5867 

0.5764 

0.5867 

0.8116 

0.0424 

0.0841 

0.6220 

0.0841 

0.8181 

0.4837 

0.6220 

0.4837 

0.6596 

0.3354 

0.7005 

0.3154 

0.7005 

0.1068 

0.7164 

0.3154 

0.7164 

0.7150 


A{:,-.,1,2) 

A{:,-.,2,1) 

A{:,-.,2,3) 


/ 0.2373 0.5028 0.7260 \ 
0.5028 0.5211 0.4278 , 

V 0.7260 0.4278 0.6791/ 

/ 0.2373 0.5028 0.7260 \ 
0.5028 0.5211 0.4278 , 

V 0.7260 0.4278 0.6791/ 

/ 0.1734 0.5972 0.6791 \ 
0.5972 0.0605 0.4080 , 

V 0.6791 0.4080 0.6569/ 
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^(:,:,3,3) = 


e(:,:,l,2) = 


e(:,:,2,l) = 


e{:,:,2,3) = 


e(:,:,3,2) = 


0.0424 0.0841 0.6220 
0.0841 0.8181 0.4837 
^0.6220 0.4837 0.6596^ 

'0.4897 0.6299 0.6104' 
0.6299 0.0527 0.5803 
^0.6104 0.5803 0.5479, 

'2.3955 2.2026 2.8921' 
2.2026 2.8852 2.5060 
^2.8921 2.5060 2.2619, 

/ 2.3955 2.2026 2.8921' 
2.2026 2.8852 2.5060 
\ 2.8921 2.5060 2.2619, 

'2.6280 2.1537 2.2689' 
2.1537 2.9841 2.2698 
^2.2689 2.2698 2.1981, 

/ 2.6280 2.1537 2.2689' 
2.1537 2.9841 2.2698 
\ 2.2689 2.2698 2.1981, 


^(:,:,3,2) = 


e(:,:,l,l) = 


e(:,:,l,3) = 


e(:,:,2,2) = 


e(:,:,3,l) = 


e(:,:,3,3) = 


0.1734 0.5972 0.6791 
0.5972 0.0605 0.4080 
^0.6791 0.4080 0.6569, 

/ 2.5328 2.6133 2.7630' 
2.6133 2.5502 2.4151 
\ 2.7630 2.4151 2.3012, 

'2.5186 2.8867 2.7372' 
2.8867 2.4538 2.2579 
^2.7372 2.2579 2.1332, 

'2.9037 2.7948 2.5391' 
2.7948 2.1978 2.2653 
^2.5391 2.2653 2.4799, 

'2.5186 2.8867 2.7372' 
2.8867 2.4538 2.2579 
^2.7372 2.2579 2.1332, 

'2.9427 2.3596 2.7611' 
2.3596 2.7011 2.6822 
,2.7611 2.6822 2.6665, 


Note that the stopping criterion in Algorithm 1 is || = 0 for exactly solving 
TGEiCP. In practical implementation, we usually use 


RelErr : = := || < Tol 


(5.4) 


as the termination criterion to pursue an approximate solution with a preset toler¬ 
ance ‘Tol’. Now, we test three scenarios of ‘Tol’ by setting Tol := {5 • 10“^, 10“^, 
5 • 10“^}. We consider two cases of the starting point where the first case is 
a vector of ones, i.e., = (1, • • • , 1)^, and the second one is a random vector 

uniformly distributed in (0,1), (the corresponding Matlab script is randCn,!)). 
To demonstrate the reliability of Algorithm 1, we report the number of iterations 
(‘Iter.’), computing time in seconds (‘Time’), the relative error (‘RelErr’) dehned 
by (5.4), eigenvalue (‘EigValue’) and the corresponding eigenvector (‘EigVector’). 
The computational results with respect to different initial points are summarized 
in Tables 1 and 2, respectively. 

From the data reported in Tables 1 and 2, it is clear that our Algorithm 1 can 
successfully solve the TGEiCP, even though it seems that the number of iterations 
increases significantly as the decrease of tolerance ‘Tol’. Actually, we tested a series 
of random starting points, and observed that random starting points often perform 
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Table 1 Computational results with starting point (1, ■ ■ ■ , 1)^- 

Example 

Tol 

Iter. 

Time 

RelErr 

EigValue 

EigVector 

Example 5.1 

5.0e-03 

657 

0.17 

5.007e-03 

0.4859 

(0.2697,0.6407)^ 

Example 5.2 

5.0e-03 

231 

0.06 

5.006e-03 

1.5609 

(0.2168, 0.1532, 0.8774)T 

Example 5.3 

5.0e-03 

536 

0.14 

5.005e-03 

0.2189 

(0.0630, 0.0000, 0.7236)^ 

Example 5.1 

l.Oe-03 

3211 

0.78 

l.OOOe-03 

0.4850 

(0.2601, 0.6512)^ 

Example 5.2 

l.Oe-03 

1703 

0.44 

l.OOle-03 

1.5512 

(0.2194, 0.1576, 0.8683)^ 

Example 5.3 

l.Oe-03 

2584 

0.65 

l.OOOe-03 

0.2173 

(0.0542, 0.0000, 0.7319)^ 

Example 5.1 

5.0e-04 

6367 

1.61 

5.001e-04 

0.4849 

(0.2589,0.6525)^ 

Example 5.2 

5.0e-04 

2929 

0.77 

5.002e-04 

1.5514 

(0.2199, 0.1575, 0.8678)^ 

Example 5.3 

5.0e-04 

5293 

1.52 

5.000e-04 

0.2171 

(0.0530, 0.0000, 0.7330)^ 


Table 2 Computational results with a random starting point. 


Example 

Tol 

Iter. 

Time 

RelErr 

EigValue 

EigVector 

Example 5.1 

5.0e-03 

277 

0.09 

5.005e-03 

0.4859 

(0.2697, 0.6407)^ 

Example 5.2 

5.0e-03 

291 

0.09 

5.003e-03 

1.5464 

(0.2172, 0.1600, 0.8673)T 

Example 5.3 

5.0e-03 

519 

0.14 

5.003e-03 

0.2189 

(0.0623, 0.0008, 0.7234)T 

Example 5.1 

l.Oe-03 

3218 

0.80 

l.OOOe-03 

0.4850 

(0.2601, 0.6512)T 

Example 5.2 

l.Oe-03 

1613 

0.43 

l.OOle-03 

1.5511 

(0.2195, 0.1577, 0.8680)^ 

Example 5.3 

l.Oe-03 

2636 

0.70 

l.OOOe-03 

0.2173 

(0.0542,0.0000, 0.7319)^ 

Example 5.1 

5.0e-04 

6071 

1.54 

5.000e-04 

0.4847 

(0.2565,0.6551)^ 

Example 5.2 

5.0e-04 

2341 

0.64 

5.002e-04 

1.5510 

(0.2203, 0.1576, 0.8672)T“ 

Example 5.3 

5.0e-04 

5341 

1.41 

5.001e-04 

0.2171 

(0.0530,0.0000,0.7330)^ 


better than the deterministic vector of ones in terms of taking less iterations as 
reported in Table 2. However, all experiments show that Algorithm 1 is reliable 
for solving TGEiCP. 

Taking a revisit on Algorithm 1, the iterative scheme (5.2) plays an signifi¬ 
cant role in the whole algorithm. In other words, the projection step given in (5.2) 
dominates the main task of Algorithm 1. As we know, the typical projection meth¬ 
ods consist of two important components, i.e., step size and search direction. In 
Algorithm 1, and serve as the step size and search direction, respectively. 
It is well known that good choices of step size and search direction may lead to 
promising numerical performance. Turn our attention to (5.2), it can be easily 
seen that step size approaches to zero as the sequence gets close to a 

solution of TGEiCP, thereby reducing the speed of convergence of Algorithm 1. A 
naturally simple idea is to increase Sk by attaching a larger constant a to it, that 
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is, the projection step in (5.2) turns out to be 

^(fc) _ ]j^ , (5,5) 

In our experiments, we observe that Algorithm 1 could be accelerated greatly when 
we set a £ (1,8). We also report some computational results in Table 3. 


Table 3 Computational results with starting point (1, • • • , 1)^ and o = 5 in (5.5). 


Example 

Tol 

Iter. 

Time 

RelErr 

EigValue 

EigVector 

Example 5.1 

5.0e-03 

130 

0.04 

5.012e-03 

0.4859 

(0.2696,0.6407)^ 

Example 5.2 

5.0e-03 

62 

0.02 

5.006e-03 

1.5472 

(0.2168, 0.1597, 0.8682)T 

Example 5.3 

5.0e-03 

105 

0.03 

5.010e-03 

0.2189 

(0.0629, 0.0000, 0.7236)^ 

Example 5.1 

l.Oe-03 

639 

0.17 

l.OOle-03 

0.4850 

(0.2601,0.6512)T 

Example 5.2 

l.Oe-03 

230 

0.07 

l.OOle-03 

1.5501 

(0.2204, 0.1580, 0.8664)T 

Example 5.3 

l.Oe-03 

513 

0.13 

1.001e-03 

0.2173 

(0.0542, 0.0000, 0.7319)T 

Example 5.1 

5.0e-04 

1270 

0.32 

5.002e-04 

0.4849 

(0.2589,0.6525)^ 

Example 5.2 

5.0e-04 

549 

0.14 

5.003e-04 

1.5513 

(0.2207, 0.1574, 0.8669)^ 

Example 5.3 

5.0e-04 

1054 

0.28 

5.003e-04 

0.2171 

(0.0530,0.0000,0.7330)^ 

Example 5.1 

l.Oe-04 

6297 

1.55 

l.OOOe-04 

0.4848 

(0.2579,0.6536)^ 

Example 5.2 

l.Oe-04 

3227 

0.84 

l.OOOe-04 

1.5520 

(0.2203, 0.1571, 0.8679)^ 

Example 5.3 

l.Oe-04 

6332 

1.65 

1.000e-04 

0.2170 

(0.0518, 0.0005, 0.7337)^ 


By comparing the results in Tables 1 and 3, it is apparent that the rehned 
projection step (5.5) outperforms the original one in (5.2) in terms of taking much 
less iterations. In Fig. 1, we further consider two different projection steps, and 
graphically plot the evolutions of the relative error defined by (5.4) in the logarith¬ 
mic sense, i.e., log(||i/^^^ ||), with respect to iterations, where the stopping tolerance 
‘Tol’ is set to be Tol : = 10“'^. 

It is clear from the above results that attaching a relaxation factor a in (5.5) is 
necessary to improve the numerical performance of our algorithm. In future work, 
we will introduce a self-adaptive strategy to adjust the relaxation factor a for an 
acceleration of the proposed method. 


6 Conclusions 

This paper considers the TGEiCP with symmetric structure, which is an interest¬ 
ing generalization of matrix eigenvalue complementarity problem. To the best of 
our knowledge, the development of TGEiCP is in its infancy and such a problem 
has been received much less attention. In this paper, we discuss the existence of the 
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Fig. 1 Evolutions of ‘RelErr’ defined by (5.4) with respect to iterations. The left plot corre¬ 
sponds to the original projection scheme, i.e., a = 1. The right one is corresponding to (5.5) 
with 0 = 5. 

solution of TGEiCP under some conditions, in addition to presenting two equiva¬ 
lent optimization reformulations for the purpose of analyzing the upper bound of 
cone eigenvalues of tensors. The bounds of the number of eigenvalues of TGEiCP 
are also presented. Einally, we develop a hrst-order projection method which might 
be a better candidate for TGEiCP than second-order solvers. Note that we only 
consider the optimization reformulations of symmetric tensors, and many prob¬ 
lems are lack of such a symmetric structure. Hence, our future work will further 
study TGEiCPs in absence of symmetric property. On the other hand, our numer¬ 
ical simulations show us that the attached a in (5.5) is important for algorithmic 
acceleration. Then, how to improve the numerical performance of Algorithm 1 is 
also one of our future concerns. 
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